# constants used for the result of line and line-segment intersection tests.
PARALLEL_LINES = 0
EQUAL_LINES = 1
INTERSECTING_LINES = 2

OVERLAPPING_SEGMENTS = 3
INTERSECTING_SEGMENTS = 4

def make_bbox(x0, y0, x1, y1):
    if x0 > x1:
        x0, x1 = x1, x0
    if y0 > y1:
        y0, y1 = y1, y0
    return (x0, y0, x1, y1)

def in_bbox(x, y, x0, y0, x1, y1):
    '''
    Checks if point (x, y) is in the rectangle defined by (x0, y0) and
    (x1, y1). Points on the boundary are considered in the rectangle
    too.
    
    Returns True if the point is in the rectangle, False otherwise.
    '''
    # order the x-coordinates such that x0<=x1
    if x0>x1:
        x0, x1 = x1, x0
    # order the y-coordinates such that y0<=y1
    if y0>y1:
        y0, y1 = y1, y0
    return x0 <= x and x <= x1 and y0 <= y and y <= y1



def overlapping_bbox(bb0, bb1):
    x0, y0, x1, y1 = bb0
    x2, y2, x3, y3 = bb1
    
    return (((x0 <= x2 and x2 <= x1) or
             (x0 <= x3 and x3 <= x1) or
             (x2 <= x0 and x0 <= x3) or
             (x2 <= x1 and x1 <= x3)) and
            ((y0 <= y2 and y2 <= y1) or
             (y0 <= y3 and y3 <= y1) or
             (y2 <= y0 and y0 <= y3) or
             (y2 <= y1 and y1 <= y3)))

def check_intersection(x0, y0, x1, y1, x2, y2, x3, y3):
    '''
    Checks if the lines through (x0,y0)-(x1,y1) and (x2,y2)-(x3,y3) intersect.
    Based on http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/

    The coordinates x0--y3 shall be floats (integer coordinates are more likely
    to result in rounding errors).

    Returns a pair (result, point), where result is one of
    the constants defined above:
    
      PARALLEL_LINES
      EQUAL_LINES
      INTERSECTING_LINES
      OVERLAPPING_SEGMENTS
      INTERSECTING_SEGMENTS
      
    If result is INTERSECTING_LINES or INTERSECTING_SEGMENTS, the coordinates
    of the intersection can be computed with get_intersection()
    '''
    result = PARALLEL_LINES
    
    x0_x2 = x0 - x2
    x1_x0 = x1 - x0
    x3_x2 = x3 - x2
    
    y0_y2 = y0 - y2
    y1_y0 = y1 - y0
    y3_y2 = y3 - y2
    
    denom = (y3_y2 * x1_x0) - (x3_x2 * y1_y0)

    noma = ((x3_x2 * y0_y2) - (y3_y2 * x0_x2))
    nomb = ((x1_x0 * y0_y2) - (y1_y0 * x0_x2))
    
    if denom == 0:
        if noma == 0 and nomb == 0:
            result = EQUAL_LINES
            bb0 = make_bbox(x0, y0, x1, y1)
            bb1 = make_bbox(x2, y2, x3, y3)
            if overlapping_bbox(bb0, bb1):
                result = OVERLAPPING_SEGMENTS
        else:
            result = PARALLEL_LINES
    else:
        result = INTERSECTING_LINES
        ua = noma / denom
        ub = nomb / denom
        if (0 < ua) and (0 < ub) and (ua < 1.0) and (ub < 1.0):
            result = INTERSECTING_SEGMENTS

    return result

def get_intersection(x0, y0, x1, y1, x2, y2, x3, y3):
    '''
    Compute the intersection of the line through (x0,y0)-(x1, y1)
    and the line through (x2,y2)-(x3,y3).
    Assumption: the lines have exactly 1 intersection point.
    If the assumption does not hold, the result is undefined.
    
    Use check_intersection to decide if the assumption holds:
    it holds when check_intersection returns INTERSECTING_LINES
    or INTERSECTION_SEGMENTS.
    '''
    try:
        x0_x2 = x0 - x2
        x1_x0 = x1 - x0
        x3_x2 = x3 - x2
        
        y0_y2 = y0 - y2
        y1_y0 = y1 - y0
        y3_y2 = y3 - y2
        
        denom = (y3_y2 * x1_x0) - (x3_x2 * y1_y0)
    
        noma = ((x3_x2 * y0_y2) - (y3_y2 * x0_x2))
        
        ua = noma / denom
        
        xa = x0 + ua*(x1_x0)
        ya = y0 + ua*(y1_y0)
      
    except ZeroDivisionError:
        xa = ya = None
    return (xa, ya)

def overlapping_polygons(poly0, poly1):
    '''
    Checks if two polygons overlap. poly0 and poly1 are lists
    of 2-tuples of coordinates, e.g., [(x0,y0),(x1,y1),(x2,y2)].
    Coordinates can be given as integer or float.

    Polygons overlap if there is at least one point that lies on at
    least one edge of each polygon.
    
    Returns True if the polygons overlap, False otherwise.
    '''

    for x,y in poly0:
        if inside_polygon((x,y), poly1, edge=True):
            return True

    for x,y in poly1:
        if inside_polygon((x,y), poly0, edge=True):
            return True

    return False

def poly_bbox(poly):
    x0, y0 = x1, y1 = poly[0]
    for x,y in poly:
        if x<x0:
            x0 = x
        if x>x1:
            x1 = x
        if y<y0:
            y0 = y
        if y>y1:
            y1 = y
    return (x0,y0,x1,y1)

def count_intersections(line, poly):
    '''
    Checks if point p is inside polygon poly. If p is on a line segment of poly,
    it is considered inside poly as well.
    Assumption: all coordinates are non-negative.
    '''
    x0, y0 = line[0]
    x1, y1 = line[1]
    x2, y2 = poly[0]
    # add closing segment to poly
    cpoly = poly + [(x2, y2)]
    
    inside = 0
    for x3, y3 in cpoly[1:]:
        v = check_intersection(x0,y0,x1,y1,x2,y2,x3,y3)
        if v == INTERSECTING_SEGMENTS:
            inside += 1
        elif v == OVERLAPPING_SEGMENTS:
            inside = len(cpoly)
            if inside % 2 == 0:
                inside += 1
            return inside 
        x2, y2 = x3, y3
    return inside

def in_closed_interval(x0, x, x1):
    return x0<=x and x<=x1

def point_on_line(x, y, x0, y0, x1, y1):
    return (x1-x0)*(y-y0) == (y1-y0)*(x-x0)
            
def point_on_segment(x, y, x0, y0, x1, y1):
    return ((in_closed_interval(x0, x, x1) or in_closed_interval(x1, x, x0))
            and (in_closed_interval(y0, y, y1) or in_closed_interval(y1, y, y0))
            and point_on_line(x, y, x0, y0, x1, y1)
            )

def inside_polygon(p, poly, edge=False):
    '''
    Computes if a point is inside a polygon.
    Copied from C-code on: http://www.visibone.com/inside_polygon/
    Adapted with the "edge" flag to handle points on polygon edges.
    '''
    inside = False
    npoints = len(poly)
    if npoints < 3:
        return False;
    xt, yt = p
    xold=poly[-1][0]
    yold=poly[-1][1]
    for xnew, ynew in poly:
        if point_on_segment(xt,yt,xold,yold,xnew,ynew):
            return edge
        if (xnew > xold):
            x1 = xold
            x2 = xnew
            y1 = yold
            y2 = ynew
        else:
            x1 = xnew
            x2 = xold
            y1 = ynew
            y2 = yold
        if  ((x1 <= xt)
            and (xt < x2) # edge "open" at one end
            and (yt-y1)*(x2-x1) < (y2-y1)*(xt-x1)):
            inside = not inside
        xold=xnew
        yold=ynew
    return inside


def point_in_polygon(p, poly):
    '''
    Checks if a point is inside a polygon (or on an edge or vertice of the 
    polygon)
    This function used count_intersections() to compute the number of 
    intersection points of a horizontal line-segment starting at p going infinitely
    far to the right (positive x-axis) and all line segments of the polygon.
    
    The function fails if the horizontal line-segment passes through polygon vertices. 
    '''
    pbbox = poly_bbox(poly)
    x0, y0 = p
    # as a possible optimization, we could check if p is outside pbbox
    
    # segment (x0,y0)-(x1,y1) has at least 1 endpoint outside poly
    x1, y1 = (pbbox[2] + 1.0, y0)
    line = [(x0,y0), (x1,y1)]
    intersections = count_intersections(line, poly)
    return intersections % 2 == 1

def winding_number(p, poly):
    '''
    Computes the "Winding Number" for polygon p about point p.
    '''
    if len(poly)<3:
        return 0
    
    # translate poly so p is the origin
    transp = [ (x-p[0], y-p[1]) for (x,y) in poly ]
    wnum = 0.0
    x0, y0 = transp[-1]
    for x1, y1 in transp:
        if y0*y1 < 0.0: # segment (x0,y0) - (x1, y1) crosses the x-axis
            
            # compute x-coordinate of intersection with x-axis
            r = x0 + ((y0 * (x1-x0))/(y0-y1))
            
            if r > 0.0: # intersects at positive x-axis
                if y0 < 0.0: # in upwards direction
                    wnum += 1.0
                else: # in downward direction
                    wnum -= 1.0
        elif (y0 == 0.0) and (x0 > 0.0): # (x0, y0) on positive x-axis
            if (y1 > 0.0):
                wnum += 0.5
            else:
                wnum -= 0.5
        elif (y1 == 0.0) and (x1 > 0.0): # (x1, y1) on positive x-axis
            if (y0 < 0.0):
                wnum += 0.5
            else:
                wnum -= 0.5
        # else: no change to wnum
        x0, y0 = x1, y1
    return wnum

def point_in_polygon2(p, poly, edge=False):
    x0, y0 = poly[-1]
    for x1, y1 in poly:
        if point_on_segment(p[0], p[1], x0, y0, x1, y1):
            return edge
        x0, y0 = x1, y1
    return winding_number(p, poly) != 0.0

def bbox_width(bbox):
    return abs(bbox[2] - bbox[0])

def bbox_height(bbox):
    return abs(bbox[3] - bbox[1])

def bbox_area(bbox):
    return bbox_width(bbox)*bbox_height(bbox)

def bbox_center(bbox):
    return (bbox[0] + bbox_width(bbox)/2,
            bbox[1] + bbox_height(bbox)/2)

    
if __name__ == "__main__":
    x0, y0 = 0.0, 0.0
    x1, y1 = 10.0, 10.0
    x2, y2 = 6.0, 0.0
    x3, y3 = 6.0, 6.1
    
    assert check_intersection(x0, y0, x1, y1, x2, y2, x3, y3) == INTERSECTING_SEGMENTS
    
    poly = [(-3.0, -3.0), (3.0, -3.0), (3.0, 3.0), (-3.0, 3.0)]
    w = winding_number((0.0, 0.0), poly)
    print "Winding num:", w, "(should be 1.0)" 
    
    print point_in_polygon2((0.0, 0.0), poly)
